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1 Introduction 


In recent years, interest in long-time integration of partial differential equations has increased 
greatly due to the need to solve diverse problems occurring in various fields of science and engi- 
neering such as fluid mechanics, aeroacoustics, electromagnetism and others. In such computations 
a large system of equations has to be evaluated (explicit schemes) or solved (implicit schemes) 
for many time steps. These simulations require huge computation time and unless more effi- 
cient computational methods are developed, they will be practically intractable in the foreseeable 
future. 

The two possible approaches to employ finite difference approximations to such problems 
are either to use a highly accurate scheme (i.e., a high order scheme or a scheme for long-time 
integration [10, 13]) on a grid which resolves all the physical frequencies occurring in the problem, 
or to use a low order scheme on a significantly finer grid, or a combination of these two. The 
first approach seems more appealing theoretically. Indeed, high order spatial discretizations [6, 
7, 13] and discretizations for long-time integration [10, 13], as well as high order time marching 
schemes [8] have been in the focus of research lately. Although significant progress has been made, 
the two main problems investigated in the abovementioned research still lack general solutions. 
The appropriate treatment of the boundary terms in high order Runge-Kutta schemes that will 
maintain the interior discretization accuracy still requires further investigation even for linear 
variable coefficient equations. This problem is mainly of a theoretical interest as it has only a 
minor effect on most practical computations [8]. The second problem is the lack of a systematic 
method for constructing numerical boundary conditions of a required accuracy such that the 
resulting discretization is time-stable. This is a major obstacle to long-time simulations. The 
Large Discretization Step (LDS) methods, presented here, offer a new and interesting approach 
to long-time integration. They enable to obtain a fine grid accuracy by time stepping mainly 
on a coarse grid with rare visits to the fine grid, at a cost substantially smaller than fine grid 
simulation. 

In some cases, the huge computational cost of fine grid simulations may be reduced by using 
such a grid only at regions where it is required, e.g., to resolve shocks, and employing coarser grids 
in parts of the computational domain where the solution is relatively smooth [1, 2, 3]. This method 
of local mesh refinement for systems of conservation laws has been reported to achieve a speedup 
of up to a factor of 55 for three dimensional problems, relative to performing the computation 
on a uniform grid with the finest mesh employed [1], This approach assumes the scheme has 
the same spatial and temporal accuracy and does not seem applicable to implicit schemes. The 
programming effort involved in generating and moving the fine grid patches is probably the cause 
for the limited use of this method. 

Multigrid methods have been employed to accelerate time dependent computations in several 
ways. The naive approach is to use an efficient multigrid solver for implicit time marching schemes. 
However, in this approach one is still confined to the fine grid time step. A more advanced idea 
is to use multigrid in time, as well. This approach, the frozen r method, was successfully applied 
to parabolic equations [4, 5, 9]. There, correction terms are added to the coarse grid equations 
such that one can time-step on the coarse grid and practically obtain the fine grid solution. This 
method exploits the smoothness of the change in the solution, typical to parabolic equations, 
which can be well approximated on coarser grids. 
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The Large Discretization Step (LDS) methods, first introduced by the authors in [11], may be 
viewed as a generalization of the frozen r method aimed at accelerating the solution of hyperbolic 
as well as parabolic equations, for both implicit and explicit time marching schemes. The present 
work investigates the LDS approach for hyperbolic equations; a class of equations that was not 
previously amenable to multigrid methods. Although the error bounds derived in this work apply 
to any time dependent equation, the algorithm for parabolic equations would significantly differ 
from the hyperbolic solver. Nevertheless, it is expected that introducing the ideas outlined in this 
work to parabolic solvers will substantially improve their performances, as well. 

The present paper extends the preliminary results presented in [11] in several important as- 
pects. The algorithm for high degree LDS was significantly improved; resulting in a more efficient 
algorithm that requires lower order intergrid transfers. The method was extended to treat non- 
linear problems with the same efficiency. Last, problems with non-periodic boundary conditions 
were solved. 

The hyperbolic LDS identifies two grids, a coarse representation grid on which all wavelengths 
occurring in the physical problem can be well resolved, and a finer computational grid which is 
required to obtain the desired accuracy with the given discretization at the prescribed final time. 
The LDS method employs a grid possibly finer than the representation grid, yet significantly 
coarser than the computational grid. It introduces one or more correction terms to the coarse 
grid equations and a system of equations satisfied by these terms is derived, initialized using 
the fine grid and solved on the coarse grid to yield the fine grid solution. The correction terms 
are integrated on the coarse grid, hence, their accuracy deteriorates at a rate determined by the 
coarse grid discretization. However, since their norm is significantly smaller than the solution 
norm, they can be effectively used for many coarse grid time steps. Thereafter, the fine grid 
should be revisited to compute new initial data for the correction equations. 

The LDS method assumes that a grid which resolves all the physical frequencies occurring 
in the problem as well as a discretization suitable for a fairly long simulation time are given. 
However, the requirement to employ the same discretization for substantially longer integration 
time while maintaining a desired accuracy necessitates the use of significantly finer grids. The 
algorithm solves on the coarse grid an extended system of equations, using essentially the original 
time marching subroutines (with at most slight modifications), yielding the fine grid solution. This 
programming simplicity renders the proposed method easily applicable to any problem similar to 
these investigated in this work, provided it obeys a few programming conventions. This is an 
important feature of the proposed method. 

The efficiency of the LDS is defined as the cost of computing the solution on the fine grid 
relative to the cost of obtaining the same solution with the LDS on the coarse grid. The typical 
efficiency achieved in this work was 16 for 2D problems and 28 for 3D problems. This efficiency 
was obtained for linear problems with periodic and Dirichlet boundary conditions and for the 
nonlinear Euler equation in a periodic domain. A particularly good discretization yielded, for a 
linear problem, an efficiency of 25 in 2D and 66 in 3D. 

The organization of this paper is as follows. Section 2 contains a heuristic derivation of the 
method for both linear and nonlinear time dependent equations. Section 3 presents bounds on 
the error in the LDS approximation for linear equations. In Section 4 it is shown that the LDS 
approximation maintains the stability and consistency of the original scheme. The LDS algorithms 
are described in Section 5. In Section 6 Fourier analysis is employed to analyze various aspects of 
the algorithm, in particular, to obtain the necessary orders of the intergrid transfers. Section 6 
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presents numerical results and the conclusions are given in Section 7. 


2 Heuristic Derivation 


An intuitive and informal derivation of the LDS method will be outlined in this section. A rigorous 
derivation for linear problems is given in Section 3. 

Consider a linear time dependent system of equations with coefficients possibly dependent on 
x but not on 


u t = 

L(x,£)u 

for x 6 ft 

Mu = 

0 

for x € dfl 

u(x,0) - 

u o(x) 

for x £ Q. 

d_ _ (_d d_ 

dx v dx\ ’ dx2 ’ 

' dxd / * 



*e[o,r] 


( 2 . 1 ) 


-L V / O 

Assume that the fine grid is required to obtained the desired accuracy at time T. However, instead 
of solving on the fine grid 


Ut = L h U h (2.2) 

one would like to modify the coarse grid equation such that it will yield the fine grid solution. 
The coarse grid solution satisfies, 


U t H = L h U h (2.3) 

A correction term r to equation (2.3) is sought such that, 

Ut = L H u h + r (2.4) 

where u h - I^U h denotes a restriction of the fine grid solution to the coarse grid. The relative 
error, u h - U H , satisfies 


(u h -U H ) t = L H (u h -U H ) + r (2.5) 

Thus, this error satisfies the same equation as u h . Moreover, 

(d t - L h )t = (d t - L H fu h = n (2.6) 

If the following relation holds, which is reasonable to assume when L H well approximates L h , 

(d t - L H fu k < (d t - L H )u h (2.7) 

i.e., n < r, then t x may be neglected; otherwise, the same argument implies that satisfies a 
similar equation to u h , resulting in a larger system of correcting equations (See Section 3). Thus, 
when relation (2.7) holds and r is properly initialized the system of equations 

uj 1 = L H u h + r (2.8) 

_ rtf 

T t = L T 
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yields the fine grid solution on the coarse grid for some integration time. This method of ex- 
panding a system of difference equations to obtain a more accurate approximation will be called 
an LDS approximation. In particular, an inflated system of the type (2.8) will be called an LDS 
approximation in Correction Scheme form. 

Introducing the new variable v h = u h + r, the system (2.8) can be written as, 

= L H u h + v h -u h (2.9) 

v? = L H v h + v h -u h 

This later form might look awkward; however, as will be shown next, this is the form of the LDS 
for nonlinear problems. This representation of the LDS will be called Full Approximation Scheme , 
following the multigrid naming conventions [4]. 

In Section 3 rigorous error bounds on such approximations for linear evolution equations are 
derived. 

For nonlinear problems, only a heuristic derivation is given. Consider the nonlinear evolution 
problem, 

U t = P(x,u,-^~) for x € D t € [0,T] 

OX 

M(u) = 0 for x € dQ (2.10) 

u(x,0) = Uo(x) for 

where D C M d , and £ = (577, 5^, fy- 

Let P h ,P H be the same semi-discretization of equation (2.10) on grids h and H , respectively. 
A modification of the coarse grid equation is sought that will yield on that grid the fine grid 
solution, i.e., a forcing term is required which will satisfy 

Uf = P H {u h ) + T (2.11) 

where u h = I?U h denotes a restriction of the fine grid solution to the coarse grid. In this case 
the relative error satisfies, 

( u h -U H ) t = P H (u h )- P H (U H ) + r (2.12) 

» P“{u h )(u h - U H ) + r 

where P*? (u h ) is a linearization of P H around u h . Assume that the following relation holds, which 
is reasonable if P H well approximates P h , 

(d, - P”(u h ))r = (d, - />"(«*)) V - v") < (9, - if («‘))(«* - V H ) (2.13) 

then the right hand side of the r equation may be neglected. 

It follows, by the same argument as in the linear case, that the system 

u h t = P H {u h ) + r (2.14) 

T t = P?{u k )T 

yields on the coarse grid the fine grid solution, for some integration time. 
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This formulation of the LDS is inconvenient to use as it requires explicit linearization. Intro- 
ducing the new variable = u 11 + r results in the following system of equations, 

Ut = P H (u h ) + v h - u h (2.15) 

vf = P H (v h ) + v h - u h 

In this setting, the implementation of the nonlinear LDS necessitates only minor modifications to 
the original time marching program. 


3 Approximation Theorems 

The LDS approximation introduced in Section 2 may be better appreciated once the initial con- 
ditions for the correction equations are determined, and error bounds on these approximations 
are derived. These two issues are the subject of the present section. 

First, an error bound is obtained for a semi-discrete approximation in a restricted setting. This 
restriction, a commutativity assumption, is introduced merely to maintain a simple and intuitive 
presentation. Subsequently, error bounds for semi-discrete and fully discrete approximations are 
derived without this superfluous assumption. 

The analysis will be performed for linear equations with coefficients which may depend on x 
but not on f, of the form, 


Ut = L(x,—)u + F(x) 

for x € D t € [0, T ] 


Mu = 0 

for x 6 dtt 

(3.1) 

u(a:,0) — uq(x ) 

for x € D 



where SI C and & = (£, jfj, • • • , £). 


3.1 Semi-Discrete Analysis 

3.1.1 Motivation 


Consider a stable semi-discretization of a linear homogeneous initial 
the form (3.1) given by, 


du h 

dt 


- L h u h = 0 


boundary value problem of 


(3.2) 


with initial conditions u A (0) = Uq. 

Let L h be an approximation to L h , e.g., a coarse grid representation of the fine grid operator. 
Define the system 


/ 


"0,771 


\ 


_d 

dt 


{ v T’ m ) 


with initial data 


/ L h I 0 
0 L h I 


L h ) 


\ 



( o > 



=: 

0 

/ 

\ ^m,m / 


w 


^m(O) = (l h - l h y U h (0) 


(3.3) 


(3.4) 
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Henceforth, an approximation of a system of equations by an enlarged system of the form (3.3) 
will be called an LDS approximation of degree m. 

Assume that L h and L h commute. The solution of this system of linear ordinary differential 
equations is e Aht (vft m (0 ), . . ., m (0)) r , where A h denotes the above system. The matrix A h has 

a block Jordan form; hence, an explicit expression for v^ m (t) is given by, 


m trh f h\kuk -. / . (T h — f h \ m + 

4,<o = ^ 1 E ( } “'•w = ‘ 1 '- L 11 - ( 1 { 


k = 0 


(m + 1)! 


**(o) 


for some £ G [0, /]. 

Assume that ||e^ i u /l (0)|| < C /l e /3h< ||u /l (0)||, for constants C h ,f3 h -, then 
llAO - 4.(011 = (3.5) 


The bound (3.5) implies that for any fixed final time T, limm-n*, ||w /l (T) - VQ m (T)\\ = 0, with 
convergence rate depending on the magnitude of the relative truncation error, \\L h — L h ||. 

This bound also suggests that when the relative truncation error is small, then there exists a 
time To such that the error in the approximation of a fixed degree m is small for T < To. This 
observation motivated the LDS algorithm and enables its high efficiency. In this algorithm L h 
and L h stand for the same discretization on the fine and coarse grids, respectively. The algorithm 
computes initial conditions for the correction equations using the fine grid, then marches with the 
enlarged system on the coarse grid as long as the LDS error relative to the fine grid solution is of 
the same magnitude as the error in the later solution. Then, the fine grid is revisited to compute 
new initial data for these equations. In this manner the fine grid accuracy can be obtained when 
time marching mainly on the coarse grid. 

The identity in (3.5) can be used to obtain the following inequality, 


ll« fc (0-<m(0ll< 


e £/,f || IKX' 1 - Z /l ) m+1 u ft (0)|| f m+1 

(m + 1)! 


(3.6) 


The stability of the semi-discretizations considered implies that there exists a mesh size ho and 
constants C, j3 such that for all grids with mesh size h < ho , 


||e L ' 1 '^(0)||<Ce' 3 '||Ao)|| 


(3.7) 


Therefore, for h < ho the inequality (3.6) implies that, 


iao- 4,(011 < ii( ^- z ( t; i ; ) : (o)iir+ ' ^ 


(3.8) 


Hence, for a fixed degree m and fixed integration time T, the error in the LDS approximation 
satisfies, || u h (t) - VQ m (t)\\ = 0(h^ m ^ v ), provided ( L h - L h ) = 0(h p ). Furthermore, if p > 0, 
which is reasonable to assume, this bound suggests that the LDS error decreases as M 771 * 1 )? as 
mesh is refined. 

In the next section, an error bound is proved for non homogeneous equations without the 
commutativity assumption 


6 



3.1.2 Error Bound 


Let L be a discretization of the spatial operator L in (3.1). Let A ^ be an approximation to 
and denote B* 1 = L ^ — A h . Intuitively, L ^ may be viewed as a fine grid discretization and A h 
a coarse grid approximation to L h . However, it should be noted that all the operators L h ,A h ,B h 
are defined on the same grid. The semi-discretization of equation (3.1) may be written as 

u h _ A h u h _ B h u h + fh 

u h (0) = ug ( 3 - 9 ) 

Denote the solution of this problem by U h (t) £ R N , where N is the number of grid points. 
Assume the following inequalities hold, 




< 




11^(011 < Ce^(||ttg|| + ||/ k ||) 


(3.10) 

(3.11) 

for constants C, (3 which may depend on h. The stability of the semi-discretizations considered 
implies that for fine enough grids (3.10)-(3.11) may be bounded independently of h. 

Define the system of equations 


/ 


d_ 

dt 


Vj} 

v Q,m 


\ 


v h , 

r m — 1 ,m 

yh 

m,m 


( A h 

o 


i 

A h 


0 

1 


A h 


\ 

yh \ 

^0 ,771 1 


f 0 h \ 



— 

ft - 1 

/ 

V m,m ^ 


fm / 


(3.12) 


Vtfm(O) ^ 


b£u£ \ 


= 

Bm-lUo 

v£, m (0) 


{ bM / 


(3.13) 


f- = Bjf h 

with the Bj defined inductively by 


0 < j < m 


(3.14) 


Bq = I 

B$ +1 = [ B j 1 , A h ] + B^B h j > 0 (3.15) 

It will be shown that the first component of the solution of this system, V£ m (t), approximates 

U h {t). More precisely, a bound on ||V 0 * m («) - U h {t)\\ will be presented which, for fixed t, tends to 
zero as m -+ oo. 

The vector (B^U k (t ), . . . , B^U h {t)) satisfies the equation 



BoU h \ 


A h I 0 ... \ 


( B*U h \ 


f 0 h \ 

d 

; 


0 A h I ... 


: 



dt 

B h m -xU h 




iV h 

— 

ft - 1 


B h m U h j 


* h J 


BlV h J 


fi + B h m + x U h / 
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(3.17) 


B%U h ( 0) \ 


\ 

B h m ^U\ 0) 


Bi-^Z 

B h m U h ( 0) ) 


bM / 


It will be shown that V* m (t) converges to B^U h {t) as m -*• oo. In particular, V£ m (t) converges 
to U h (t) as to — * oo. 

Define e£ m (t) = B$U h (t) - It satisfies the following equation 


_d 

dt 


( eg, m \ 


e h 

771 — 1 ,771 

V ) 


o # / ... 

A h 


\ 

/ ^ \ 


f o \ 


m — 1 ,77i 

= 

0 

/ 

\ ^m,m / 


B h m+ iU h / 


/ 4m( 0) >| 


( 0 \ 

O 

7 

= 

0 

V e m, m (0) ) 


\ 0 / 


The solution of this system satisfies 

<„(«) = 

Jo 

= f c^(*-*)ej +1>m (a)d* 0 < j < m - 1 

v 0 

The norm ||e§ m (f)|| is the sought error bound. 

Il c m,m(*)ll < l|«t+lH jf ««*— ) ||V k (*)||<fa < ||Bi + ,|| (ll u oll + ll/ ft |l) 
By induction one obtains 


(3.18) 


(3.19) 


(3.20) 

( 3 . 21 ) 


(3.22) 


f m-l + 1 

ll4.MII S n«i+ill (ll«Sll + ll/' , ll) c (m _| +1) / ‘ < 3 - 23 ) 

The following theorem was proved. 

Theorem 1 Let U h (t) be the solution of (3.9) and (V 0 * m (i), . . ., V^ m (t)) be that of (3.12)- 
(3.13)' Then 

fm + l 

itvi„(o - vHm < iiiw.ii (irtii + nrt) c ( 3 -24) 
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A different bound on the error in the LDS approximation will be derived next. Denote 

V h (t) = sup e Ps (3.25) 

0 <s<t 

*m+i(0 = sup ||i& +1 tf fc (*)|| (3.26) 

0 < 5 << 

The exponent 3 may depend on the mesh size h. However, for a grid fine enough the stability of 
the semi-discretization implies that the solution can be bounded independently of h. Therefore, 
it will be assumed that the grid is fine enough for this property to hold. Then, 

l|em,m(OII < (3.27) 

and 

r 1771 7 + 1 f 771 — 

ll e i.m(OII < [* (0] + (3.28) 

Therefore, 

ll«o,m(*)ll < [»(<)] $ "+i( f ) ( m+1 ); (3-29) 

This bound is very crude and can be easily improved. However, it implies that for a fixed final 
time T and degree m, if 5 m+1 = 0(/? Pm+1 ) with p m+1 > 0, then $* +1 (r) = 0(h p ™+') and the 
error in the LDS approximation decreases like h Pm + l as mesh is refined. 

The bound (3.29), suggests that looking at i?^ +1 may give an insight into the accuracy of the 
approximation. From the definition of Bj 1 , it follows that B ^ = B h and B% = [B h , A h ] + ( B h ) 2 . 

In the following two examples, A h is interpreted as a coarse grid approximation to the fine 
grid spatial discretization A h + B h . 


Example 1 Consider a discretization of the two dimensional wave equation 

uu = Au 

when transformed into the system 


(3.30) 



(3.31) 


This discretization of the wave equation was successfully used in [12] to solve problems in elasticity. 
Assume that the Laplace operator is approximated by a second order scheme. Then, up to higher 
order terms, 


A h — ( 0 O nh _ ( 0 0 

A + a H 2 (d XXXX + dyyyy) 0 j ’ ^ ^ — H 2 ){d XXXX -|- dyyyy ) 0 

where a denotes a generic constant. Since B 2 = 0, it follows that 

B% = [B h , A h ] — ^ — a (h ~ H 2 ) {dxxxx + dyyyy ) Q ( ^ — H 2 ) ( d XXXX + dyy yy ) \ 


(3.32) 

(3.33) 


It can be seen that B^ consists the same terms as the relative truncation error. Thus, one can not 
expect that the LDS of degree one will yield an approximation more accurate than the coarse grid 
discretization. Indeed, applying the LDS algorithm to this discretization shows that the error in 
the LDS of degree one solution grows like the coarse grid error. 


9 



Example 2 Consider the linearized Euler equation 

Pt = a ■ V + c(u x + v y ) 

Ut — cl • V + c p x (3.34) 

V t = O • V + CPy 

where a = ( 01 , 02 ) is a two dimensional vector. Assume, that the spatial operator is discretized 
using a second order scheme. Then, up to higher order terms, 


A h = 


O • (V + O. H 2 (3 rII T dyyy )) C (9x T Ot H dxxx ) 

c(d x + a H 2 d xxx ) a • (V + a H 2 (d XX x + d 

~ l a 1 u2 


C (dy + Ol H dyyy ) 

0 


B h = ( h 2 - H 2 ) 


\ 


Kyy) 

0 

a ■ (V + a H 2 

{dxxX "1" dyyy )) 

dxxx “1" dyyy ) 

C O! dxxx 

C a dyyy 


C Ot d XX x 

d * ( Ot dxxx H" dt dyyy ) 

0 

(3.35) 

C a dyyy 

0 

O * (ck dxxx dL dyyy ) J 


where a is a generic constant. It can be easily seen that ( B h ) 2 consists of sixth order mixed 
derivatives and fourth order powers of H , h. 


0 0 0 


with 


[B h , A h ] = I 0 0 T) 

0 — Tj 0 

Tj = C 2 Ol 2 (h — H ) (dyyyx — d x yyy + 9yyy x xx ~ 9 x xxyyy ) 


(3.36) 


(3.37) 

(3.38) 


For smooth solutions, 

r i =C 2 a 2 {h 2 -H 2 ) 2 (d yyy x-dxyyy ) 

Hence, the error bound (3.29) implies that for smooth data the LDS of degree one yields 
a significantly smaller error than the coarse grid operator. Indeed, our numerical results (see 
Section 7) show that the LDS algorithm for this equation yields the fine grid solution on the 
coarse grid. 

3.2 Fully Discrete Analysis 

Consider a stable finite difference approximation to (3.1) of the form 


where k = At. 
Assume that 


U n+l = (I + kA)U n + kBU n + kf 
U° = u 0 


urn < e0nk (imi + ii/ii) 

\\I + kA\\ < e 0k 


(3.39) 


(3.40) 

(3.41) 
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In an analogous way to the semi-discrete case, define (F 0 " m , . . V^ m ) to be the solution of 


/ K” +1 \ 


(I + kA kl 0 ... \ 


Km ) 


kf 0 \ 

; 


0 I + kA kl 


* 


; 

yn+l 
m — 1 ,m 




Kn n -l,m 


^ fm — 1 

V n+1 J 

\ m,m / 


\ I + kA ) 


) 


\ ^ fm J 


v° \ 

v 0,m 


f 5 0 «0 ^ 

v° , 

r m — l,77i 

= 

Bjj i — l U 0 

vim 


V B m u 0 / 

fj = BJ 

0 < j < m 


with Bj defined recursively by 

B 0 = I 

B]+ 1 = [Bj, A] 4- BjB j > 0 
Consider the vector (B 0 U n , . . . , B m - X U n , B m U n ). It satisfies 


(3.42) 


(3.43) 


(3.44) 


(3.45) 

(3.46) 


( B 0 U n+1 \ / + kA kl 0 

: 0 I + kA kl 

B m -\U n+1 
\ B m U n+1 / 


I + kA J 


( B 0 U n \ 

B m -iU n 

V B m U n j 


kf o 

k fm — 1 
V k f m + k B m+1 U n ) ) 


with the same initial condition as ( Vo m , • ..,V£ m ). 
The error e” m = BjU n - VJ satisfies 


/ e o,m \ 


n+l 
m — l,m 

\ e n+1 / 


(I + kA kl 


0 


0 

/ + £4 */ 


I + kA / 


c 0,m 


c m — 1 ,m 

e n 


V A:5 m+1 {7” ; 


(3.48) 


^ 0 ,m. 


e° 


/ 0 \ 

0 

V 0 / 


(3.49) 


(3.47) 
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These equations give, 


le” >m = '£(I + kAr- 1 -jkB m+1 U> 

3=0 


(3.50) 


+ kA) n ~ l ~i ke 3 l+l m 0 < / < m — 1 

J-0 


(3.51) 


The following theorem can be proved by bounding the solution of this discrete system. 
Theorem 2 Let U n be the solution of (3.39) and (Vo, m (0>- • • > Kn,m(<)) be that of (3-42)- 
(3-43). Then 


K m - u n \\ < (||f/°|| + 11/11) 


n U m+1 e' 9 t n_m_ i )fc 

m + 1 / 


(3.52) 


Proof : The following, more general, formula will be proved 


IKmll < (l|t , °ll + ll/ll) llftn+lll ( m _" + 1 )t"-' +1 e' , ( n - m+, - 1 >* for 0 < / < m (3.53) 
Theorem 2 is the particular case of / = 0. 

The proof follows by induction on m — l. First consider the case m - l = 0, i.e., / = m. 
According to formula (3.50) 


C,„ = EU+*- Ar-'-ikB m+1 U‘ 


(3.54) 


Therefore, 

ll e m,mll < E IIU + kAf-'-i II k IIJWi II lie'll < £ k ||Bm +1 || e Blk (|0°|l + 

j = 0 3=0 

= (||^°|| + ll/ll) ||S m +i||(nfc)^ (w - 1)fc 

Assume the bound (3.53) is correct for m - l = m - 1; that is, 

IKmll < (lie°n + ii/ii) 1^.1; ^* 

From relation (3.53) it follows that, 


eo,m = T,( I + kA ) n ~ 1 ~ jke 1 
3=0 


Therefore, by the induction hypothesis (3.56), 


KJI < £ ||(/ + fcA) n - 1 -q|* (||«7°|| + ll/ll) \\B m+1 \\( J )k m e‘ 


< (||C°|| + ll/ll) IIA 


0 (j-m) k 


n j ^ m + 1 g|0(n-m-l) k 

m+ 1) 
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The last inequality follows from the following identity which can be proved by induction on n. 

n— 1 


E 

j=0 


J 

, m , 


n 

,m + 1 


The relations (3.50)-(3.51) give rise to a somewhat different bound, as well. Denote, 

® h (n ) = sup ||(/ + &j 4) j || 


0 <j<n 


(3.60) 


(3.61) 



$m+lO) = SUp 

0<j< n 

(3.62) 

Then, 


K,«H < (n*)®V)*m+l(») 

(3.63) 

and 


Kmll < (»fc) m+1 [^(«)J m+ <+1 (n) 

(3.64) 


The stability of the discretization implies that for a fixed time T, if the meshsize is fine enough, 
i.e., h < ho , (assuming that k is related to h in a fixed manner), then one can bound 


9 h (n) < C{T ) 


for all nk < T , 


such that h < ho 


(3.65) 


this bound is independent of the mesh size. Thus, 


HeS.mll < r m+I [C(r)] m+1 < +1 (n) (3.66) 

If Bm + 1 = 0{h Prn + l ), then as mesh is refined the error in the LDS approximation decreases like 
h Pm + 1 . 

At first glance the discretization (3.39) is a first order in time explicit scheme. However, any 
single stage explicit or implicit discretization of any order can be brought to a similar form once the 
source term is appropriately redefined. Therefore, Theorem 2 can be modified and generalized 
to yield a similar result for a general discretization. Such a generalization would be futile as 
finding an explicit form for A, B and B m might require matrix inversion or computing a matrix 
polynomial. Thus, in those cases the bound (3.52) is hard to compute. 

In the fully discrete setting, as in the semi-discrete, an inspection of B m +\ might indicate 
about the applicability of the LDS method for a given discretization. 


4 Stability, consistency and convergence 

In the previous section a method for obtaining highly accurate approximations using an expanded 
system of lower order approximations was analyzed. It should be proved that if the original scheme 
was consistent and stable then so is the resulting LDS scheme. 

For simplicity of presentation the discussion is limited to LDS of degree one. The generalization 
to LDS of general degree m is straightforward. It is further assumed that the source term F = 0. 
This assumption does not effect stability analysis [14] and its effect on consistency will be shortly 
discussed. 
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In an LDS of degree one instead of using the stable and consistent scheme, 


u n+1 = ( I + kA)U n 

(4.1) 

O 

II 

es 

o 

(4.2) 

the following LDS system is employed, 


V 0 n i +1 = (I + kAW^ + kV^ 

(4.3) 

= V+* A )K 1 

(4.4) 

o 

II 

(4.5) 

I'Z = B x u o 

(4.6) 


where B x = B is defined in (3.39). 

The stability of the scheme (4.3)-(4.4) follows from the structure of this system. It consists of 
a principal part which is coupled through lower order terms. The principal part is diagonal; thus, 
its discretization is stable for the same time step as the original scheme. The lower order term 
does not affect stability. 

The consistency of the single equation discretization implies that the LDS is a consistent 
approximation of the system 

u t = Lu + t (4.7) 

T t = Lt 

for initial data 

u(x,0) = u 0 {x) (4.8) 

r(x,0) = r 0 (x) 

In the LDS approximation V^j = B x u 0 . Thus, if the following relation holds 

Bi = 0(h Pl + k qi ) with P\,q\ > 1 (4-9) 

then the LDS solution is a consistent approximation of (4.7)-(4.8) with initial data 

u(x,0) = tio(z) (4-10) 

r(x,0) = 0 

Thus, it consistently approximates the equation 

u t = Lu (4-11) 

u(x,0) = uo(z) 

The source term for the V X)1 equation is B x F. Thus, by the above reasoning, if (4.9) holds then 
the t equation is a consistent approximation of an equation with F = 0. Therefore, this term 
does not effect consistency, either. 

It follows, by Lax equivalence theorem [14], that the consistency and stability of the origi- 
nal equation together with condition (4.9) ensure the LDS solution convergence to the analytic 
solution. 
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The LDS approximation maintains the stability of the original scheme, however, it does not 
necessarily preserve the time-stability of the discretization, i.e., the LDS solution might exhibit a 
non-physical growth in time although the underlying discretization did not allow such a growth. 
This phenomena is demonstrated in Lemma 1. Recall, if A h is a discrete operator acting on 
grid functions on a periodic or infinite domain, then the symbol of A h , A h (9), is defined by the 
identity, 

A h (0)e i6x = A h e ie ' x (4.12) 

Lemma 1: Let A h be a space discretization of a scalar time-dependent equation with periodic 
boundary conditions , such that its symbol satisfies A h (0 o ) = 0, for 6 0 / 0 and B h (9 0 ) ^ 0. Then 
the error in the LDS system based on A h grows polynomially in time for this Fourier component. 
The order of the polynomial equals the LDS degree. 

Proof: Consider a semi-discrete LDS approximation of degree one based on A h . According 
to Section 3.1.1 it has the form, 

Ut = A h u h + r h (4.13) 

T t h = A h T h 

with initial data 

u h (x,0) = u£(x) (4.14) 

T h (x, 0) = B h u%(x) 

The solution of this system for the 6q component is, 

^(0o,O = ^ h($o)t (Mo(^o) + ^ /l ^o)fi ,l (^o)) = (l + tB h (6 0 )) u h (0 o ) (4.15) 

The proof for higher degree LDS is similar. 

Note, that the property A h (9o) = 0 is common to central discretizations of the first derivative 
on non-staggered grids, for 9 0 - n (see also Sec. 6.3). 

It should be emphasized that although the LDS transformation preserves stability, the resulting 
algorithm might not be stable. Consider, for instance, a discretization satisfying the condition of 
Lemma 1; then the LDS solution for the 9q component will grow exponentially with the number 
of visits to the fine grid. Other possible sources for such a growth are large errors introduced 
by the intergrid transfer employed by the algorithm which are not damped during the cycle (see 
Sec. 5.1.4) or improper use of Richardson extrapolation (see Sec. 5.2). These last remarks can be 
understood once the LDS algorithm is presented, in the next section. 

5 Large Discretization Step (LDS) Methods 

The LDS approximation, introduced in the previous sections, approximates a high accuracy 
scheme by an enlarged system of equations of a lower accuracy discretization. In the present 
work the two schemes are the same discretization of a differential operator on two different grids. 
An LDS algorithm computes the fine grid solution on the coarse grid by solving there an extended 
system of equations which are initialized using the fine grid. The LDS system is integrated on 
the coarse grid; hence, the accuracy of the correction terms deteriorates at a rate determined by 
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that grid discretization. However, since the norm of these terms is significantly smaller than the 
solution norm, they can be effectively used for many coarse grid time steps. Then, the fine grid 
has to be revisited to compute new initial data for the correction terms. These terms initialization 
is computationally costly and the LDS approximation consists of more equations than the original 
problem; however, the large number of steps performed on the coarse grid before revisiting renders 
the resulting algorithm very efficient. 

This section details the algorithmic implementation of these ideas. The algorithm is presented 
both in its Correction Scheme and Full Approximation Scheme forms and efficient initialization 
procedures for these schemes are described and analyzed. Richardson extrapolation is introduced 
to the LDS method which, for smooth solutions, yields a higher order approximation. The effi- 
ciency of the algorithm is discussed and evaluated. 

Given a system of hyperbolic differential equations with coefficients which may depend on x 
but not on f, of the form 

Mu(x, t) = 0 
U(x, 0) = U 0 (x) 

where ft C JR d , and • • •> 

Consider a discretization of the form 

U n+1 = E(x, k, h)U n + S(x, k, h)F n (5.2) 

where h,k denotes Ax and At, respectively, and U n = Uf approximates U(jh,nk). In this work 
E(x,k,h) is an explicit or implicit two level time marching operator. However, the method may 
be used with multilevel integration schemes, as well. In the sequel, the notation E k ' h will be used, 
omitting the possible dependence on x. In an LDS application two grids (in space-time) are given, 
a fine one with spacing ( h,k ) and a coarse with ( H,K ), where H = ah,Ii = ak. Given U H (x,0) 
on the ( H , K) grid, one needs to calculate the solution up to a prescribed final time T and obtain 
the fine grid accuracy. 

5.1 The LDS Method of General Degree 

The degree of an LDS approximation is defined as the number of its correction terms. The 
error bounds obtained in Section 3 suggest that in many cases (e.g., for smooth solutions) the 
higher the degree the better the LDS approximates the fine grid solution. In practice, efficiency 
considerations limit the degree to at most two (see Section 5.3.2). In the sequel, algorithms for 
LDS of general degree are described. 

The initialization procedures necessitate the transfer of the solution from the coarse grid to 
the fine and back. In this section, it is assumed that appropriate intergrid transfers are given. 
The order and properties of these transfers will be discussed in Section 6.1. For presentation 
simplicity the algorithm is described for the case F(x) = 0; the treatment of a source term is 
straightforward. 


for a: € ft, t G [0, T] 

for x G dft (5.1) 

for a; € ft 
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5.1.1 Correction Scheme LDS 


The fundamental idea of the LDS method is to look for correction terms to the coarse grid 
equations, derive and solve the equations satisfied by these terms to obtain the fine grid accuracy 
on the coarse grid. In Section 3, it was shown that for linear problems these terms approximately 
satisfy the same equation as the solution. The resulting system of equations, which is valid only 
for linear problems, was named Correction Scheme LDS. 

Correction Scheme LDS of Degree One when H = 2 ft. For clarity of presentation, the 
algorithm is first presented its simplest form, i.e., Correction Scheme LDS of degree one with 

H A_ o 

h k 

The algorithm consists two stages : initialization of the correction terms using the fine grid, 
and time marching on the coarse grid for a predetermined number of steps. The results in Section 3 
show that for linear problems the correction terms satisfy approximately the same equation as the 
solution. However, they do no indicate how to effectively and efficiently compute initial values for 
these terms. The requirement that the LDS solution should yield the fine grid solution suggests 
that on the first time steps these solutions should be identical. This observation leads to the 
following initialization routine, 

Initial ize(F 0 ^, ) 


Set 

II 

£ 

rh,k i rN 

V/C v o,i 

Solve 

u N *f = 

E h ' k U N+ ~ f 1 . 


= 

e h ' k v 0 A ; 

Set 

= 

r?f uN+1 - Vo4 +1 


& 

+ 

II 

v 0 y> + 

Set 

N 

N + 1 


The initialization consists of interpolating the LDS solution to the fine grid and stepping for 
the same time on the fine and coarse grids. Then, the correction term is set to the discrepancy 
between the two solutions and the LDS solution is updated to the fine grid. Here, Vjj A+1 stands 
for an intermediate value assigned to this variable. 

The time advance of the LDS has the following simple form, 

LDS Method of Degree One, Correction Scheme 

Initialize Vq’j 

IV = 0 

While N < [ Do 

Call Initialized^, Vj]j) 

For i =1 ,... .Revisit Do 
Solve V ^ 1 = E h > k Vg 
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v 0 a ; +i = e h - k Vqj + v^ +i 

Set N = N + 1 

End 

End 

The relation between the theory derived in Section 3 and the initialization procedure is ana- 
lyzed in Section 5.1.2. 

Correction Scheme LDS of General Degree when H = 2 ft. The correction term in the 
LDS of degree one is initialized to equalize the fine and LDS solutions. If one could solve the 
exact equation satisfied by this term, which is approximately the fine grid equation, then the fine 
grid and the LDS solutions would be identical. Clearly, this would be as difficult as obtaining 
the fine grid solution on the coarse grid. Instead, one can add a new term to correct the first 
correction term equation. Thus, in an LDS of general degree, the i th correction term may be 
viewed as correcting the (i — l) t/l term. 

The general degree algorithm consists of two stages : initialization of the correction terms 
and time marching on the coarse grid for a predetermined number of steps. The time marching 
procedure has the following simple form, where term i corrects term i — 1 for 1 < i < d. 

LDS Method of General Degree - d , Correction Scheme 

Initialize T 0 ° d 

N = 0 

While N < Do 

Call Initialized^, ...,Fj^,d ) 

For i = 1,. .. .Revisit Do 
Solve V"/ 1 = E«* V» 

For 1 = Step = -1, Do 

Solve Vfi + V™ 

End 

Set N = N + 1 

End 

End 

In an LDS of degree d, the i th correction term (1 < i < d) corrects the (i - l) </l term and both 
satisfy approximately the same equation as the coarse grid solution. Therefore, the initialization 
of the i th term to correct the (t - l) t;! term is identical to the initialization of the first term to 
correct the solution in the LDS of degree one. Once Vi t d is initialized and Vi-\4 is updated using 
this value; the lower index variables can be time advanced using these new values. This procedure 
is repeated for all 1 < i < d. Section 5.1.2 outlines the connection between the error bound 
derived in Section 3.1.1 and these intuitive arguments which are implemented in the procedure 
listed below. 

Initialized^, ...,V£ d ,d ) 
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m=l ,2 


For i = 1,. . . ,d Do 


Set 

u N = 

/M y.w , 

1 H,I< V i-l,d 

Solve 

u N+ f = 

E hM UN+Xfl' 


= 

e h.k 

Set 

v>T = 

- ift+j 


v. N+ \ = 

t — l,rf 

K-. + ] + *# +1 

For 1 


,0, Step = -1, Do 


Solve Vj +1 = E H! < V$ + 


End 

Set N = N + 1 

End 

The initialization of an LDS of degree two is graphically illustrated in Figure 1. 


U 1 = E h u 0 t 2 = E l T 1 



V-2 = Eh^\ + 


Figure 1: Initialization of an LDS of degree two. With the notation : u H = V 02 r H = V 12 

V H = v 2 , 2 
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Correction Scheme LDS of General Degree when H = ah. The above described proce- 
dure can be easily adapted for the general case when ^ = a. The simplest approach is to perform 
a time steps on the fine grid for each coarse grid time step and initialize the correction terms 
correspondingly, see Figure 2. This procedure of direct initialization is very costly and greatly 



Figure 2: Direct Initialization 

reduces the LDS efficiency. 

In case a is a composite number, e.g., a = 2 l , a more efficient approach is available, exploit- 
ing the LDS high accuracy by employing intermediate grids. In this approach, the simultaneous 
initialization , the fine grid is used to initialize an LDS system of degree one on the grid H\ = 2 h; 
since this approximation is very accurate it can be used to initialize an LDS of degree one on grid 
#2 = 2 H\. This procedure is repeated until the correction term on grid Hi = H is initialized. 
This process is repeatedly employed for all correction terms. In this method, a correction term 
on a coarse grid is initialized as soon as enough time marching on finer grids was performed, see 
Figure 3. 

A few important points should be emphasized regarding the initialization procedures. First, 
in this work the computationally efficient method was favored. There is, however, a trade off 
between the computational cost and the storage requirements of these two methods (see Sec. 5.3). 
Therefore, whenever storage is limited, direct initialization might be preferred. Second, it might 
have been expected that direct initialization will yield more accurate solutions than the simulta- 
neous approach. However, our numerical results show that these methods are indistinguishable 
for integrations times of interest, i.e., as long as the error in the LDS solution is small. Last, 
the direct initialization is of practical interest, being used with Richardson extrapolation (see 
Sec. 5.2). 

A simple way to predict the LDS performance is to look at the relative magnitude of the 
correction terms immediately after initialization. According to the result presented in Section 2, 
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Figure 3: Simultaneous Initialization 


the ratio | j |! T ^^o||| should be roughly constant. Thus, a large variation in this quantity suggests 
a large error in the initialization of Tj, causing the LDS failure. 

The error bounds derived in Section 3 apply to approximations of general degree m. In practice 
however, due to efficiency considerations two is the highest degree used (see Sec. 5.3.2). 

5.1.2 Initialization Analysis 

The procedure for the correction terms initialization was justified by the intuitive arguments that 
when properly initialized, the LDS solution should agree with the fine grid solution for the first 
time steps; and that in an approximation of degree d the i th term K\<i corrects the (i - l) th 
term i j (for 1 < i < d ). Thus, each term is initialized similarly to the manner the first 

term is set to correct the solution. This section outlines the relation between the error bound 
derived in Section 3.1.1 for the semi-discrete case with the commutativity assumption and the way 
initialization is implemented. Specifically, it will be shown that the solution at the termination of 
the initialization procedure described in Section 5.1.1 and the solution of the LDS approximation 
introduced in Section 3.1.1 are equal, up to higher order terms and multiplicative constants which 
are used for computational efficiency. 

For simplicity, the analysis is performed for the case x = 2. It is assumed that the fine 
and coarse grid spatial discretizations commute, i.e., [L H ,L h ] = 0; and that intergrid transfers 
introduce no error. By an abuse of notation these transfers are omitted from the analysis in this 
section. Nevertheless, whenever the fine and coarse grid solutions appear in the same formula, it 
should be understood that an appropriate restriction of the fine grid solution to the coarse grid 
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is implied. Similarly, whenever the fine and coarse grid operators appear in the same formula 
it stands for a restriction of the fine grid operator to the coarse grid. The properties of these 
transfers necessary to guarantee the algorithm performance are discussed .in Section 6.1. 

In Section 3.1.1 an error bound was derived for a semi-discrete LDS approximation of general 
degree d. This approximation which is given by equations (3.3)-(3.4) can be succinctly written 
as, 


dvi4 

dt 

dv^d 

dt 


= L H Vi, d + v i+ i 4 for i < 0 < d 


L H Vd,d 


with initial conditions 


v itd {x , 0) = ( L h - * u 0 (x) 

where w 0 (x) is evaluated at grid points. The solution of this system is given by, 


rd-i 






■k = 0 


e LHt uo(x ) for 0 < i < d 


(5.3) 


(5.4) 


(5.5) 


For an LDS algorithm of degree one denote by u LDS ,r the variables approximating u 0J , 
respectively. At time At the solution of the LDS approximation is given by, 


no,i («p, At) 


1 + At (L h - L H )]e LH *u 0 (x) 

e L”At ^ L h _ L H Uq(x) 


(5.6) 


For linear problems, the time marching operator E H approximates e LHAt . Therefore, in the fol- 
lowing semi-discrete analysis e LH Ai will denote the time-stepping on grid H for time A*. Assuming 
u(3, At) = 0 ( 1 ), then, at the end of the initialization phase the variables u LDS ,r satisfy 


U LDs( X i ^0 


r(x, Ai) 


e LhAt u 0 (x) = 


1 + At (i‘ - L«) + (l* - L «) 


e^ H At uo{x) + h.o.t 


v M (x, At) + 0 ((At) 2 (i‘ - r") 2 ) (5.7) 

(e 1 *" - e 1 "") «oM 

At { (l h - L”) + ^ (l* - i") 2 } e 1 " i, ««(*) + h.o.t 

At ftti,i(x. At) + O (At (z* - I") 2 ) } (5.8) 


The factor At is maintained to reduce the number of multiplications during the time-stepping 
stage. It can be seen that, up to higher order expressions, the algorithm correctly initializes the 
first correction term. 
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For an LDS algorithm of degree two denote by u lds ,t, ip the variables approximating v 0 . 2 , t>i, 2 
and v 2t 2i respectively. At time 2A t the solution of the LDS approximation is given by, 

«o, 2 (s,2A t) = l + 2At(L h -L H )+2(At) 2 (L h -L H yje LH2At u 0 (x) 
v h2 (x,2A t) = [l + 2At (L h - L h )\ e LH2At [L h - L H ) u 0 (x) (5.9) 

u 2>2 (a:,2A0 = e LH2At (z A - L H ) 2 u 0 (x) 

It can be easily seen that at the end of the initialization procedure the variables ^ LDS yT y (f satisfy. 

u LDS (x,2At) = e Lh2At u 0 (x) = v 0i2 {x,2At) + o(^At) 3 ^L h -L H y^ (5.10) 

r(x,2At) = e LhAt (e LhAt -e L * At )u 0 (x)= [ e (^-L») 2 A < _ c (L*-L»)At] e ^2At Uo(x) 

= At | ( L h -L h )+~ (j L h - X*) 2 } e LH2At u 0 {x) + h.o.t 

= At ju u (x,2Af) + 0 (^A<(L /l -L H ) 2 )| (5.11) 

<p(x,2At) = (e Lh At - e 1 ” At y u 0 {x) 

= (At) 2 1 ( L h -L h )\ At (j l h - L h ) 2 [L h + L h ) | e LH2At u 0 (x) + h.o.t. 

= (At) 2 |u 2 , 2 (a:,2At) + 0 ^At (5.12) 

It can be seen the algorithm correctly initializes the correction terms in the LDS of degree two, 
up to higher order terms and multiplicative constants. 

5.1.3 Full Approximation Scheme LDS 

The FAS form of the LDS is appropriate for both linear and nonlinear problems. However, since 
the Correction Scheme has a simpler form and necessitates less modifications to the code, it is 
more conveniently used for linear problems. 

In the present work only LDS in FAS form of degree one was implemented , and it will be 
described for a homogeneous system of equations. 

Recall, the Full Approximation Scheme of degree one is given by, 

u'l = P H (u h ) + v h - u h (5.13) 

v t h = P H (v h ) + v h - u h 

where P H may be either a linear or nonlinear operator and v h - u h 4 - r, with r corresponding to 
the first correction term in the Correction Scheme form. 

The algorithmic implementation is slightly more involved than the Correction Scheme as it 
requires some modification of the time marching procedure. Denote by E^' h the integration 
scheme obtained by modifying the coarse grid operator E H ' K to time advance the LDS system of 
degree one (5.13). 
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The algorithm has the following form 


LDS Method of Degree one , Full Approximation Scheme 

Initialize Lgj 

N = 0 

While N < rfl Do 

Call Initialize (V q^, V^) 

For i =1 , . . . , Revisit Do 
Solve V^ 1 = E?' K V t N , 

Set N = N + 1 

End 

End 

Here, V\ denotes the vector (Vo,i, Vi,i). For presentation simplicity, the initialization procedure 
is described for the case ^ = 2. Generalizations to the case ^ = 2 l are identical to those for the 
Correction Scheme. In Section 5.1.2 it was shown (see Eq. (5.8)-(5.7)), that after the Correction 
Scheme initialization Vq X contains the fine grid solution u h (x,At) and V^ t contains Atr(x, At). 
Thus, the FAS initialization consists of the Correction Scheme initialization supplemented with 
the additional computation of (u h + r)(z, At) at the end of the procedure. For completeness, the 
whole initialization procedure for the Full Approximation Scheme is listed below 

Initialized^, V^) 

Set U N = 


Solve 

II 

sl« 

+ 

to 

£ h,k jjN+'- 

71 — 1 

2 , 


yjv+l _ 
M),l - 

E h,k 

v N 

K 0,1 


Set 

II 

+ 

jH,K{jN + 1 

h,k 

1 

sH 

± 


vjT = 


+ 

«k + *. 


ViT 1 = 

v,T 

+ 

1 f/iV+l 
A t V l,l > 

Set 

N 

1 v + 

1 



The generalization to higher degree LDS is straightforward, instead of the original variables 

(Vo,<f, V\,d, • • • , V d ' d ) a new set of variables (Vo.rf, Fo,d + 1 -Vd,d) is introduced. The 

equations satisfied by these new variables can be determined from the equivalence between the 
Full Approximation Scheme and the Correction Scheme for linear problems. 

5.1.4 Treatment of Boundary conditions 

The LDS treatment of the boundary conditions will be discussed under a restrictive commutativity 
assumption, which at this stage we do not know how to dispose. 
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(5.14) 


Assume u h satisfies the boundary condition 

Mu h + Nu h = g 

If [M, A] = 0, then r = Nu h satisfies, 

Mt = M N u h = N M u h = Ng - N 2 u (5.15) 

If Mr <r the right hand side term may be neglected and r satisfies the boundary condition 



Mr = 0 

(5.16) 

Otherwise, when higher degree LDS is 

employed, the Tj satisfy the boundary conditions 



M T t +T 2 = 01 

(5.17) 


Mr k = 0 

(5.18) 

where 




t- 
' <-> 

fe; 

ii 

ir 

(5.19) 


9] = N 3 g 

(5.20) 


Thus, a large error at the boundary discretization may require adding correction terms and 
corresponding equations in the whole domain. 

The assignment of the appropriate boundary conditions to the correction terms when the 
commutativity assumption does not hold should be further investigated. 

The presence of non-periodic boundary conditions may pose problems even when the commu- 
tativity assumption holds. This is due to errors introduced by the one-sided high order interpola- 
tion near the boundaries. These large and localized errors excited during the initialization phase 
might not be damped before the next visit to the fine grid; resulting in an error which grows 
exponentially in the number of visits to the fine grid. At this stage of research, it seems that 
one should use the differential equation to design appropriate near boundary interpolation with 
reduced errors (see Sec. 7 for an example). 

5.2 Richardson Extrapolation 

The simultaneous time stepping on two grids during the initialization phase can be used to 
estimate the local truncation error. This estimate can be used in various ways. In [1, 2, 3] it was 
used to implement adaptive mesh refinement for hyperbolic equations. In the multigrid method 
this estimate is used for the r extrapolation technique which applies a weighted transfer of the 
correction term to obtain higher accuracy using a lower order scheme [4] . 

For simplicity, let Q k be a two-level explicit difference operator. If the solution is smooth 
enough, the local truncation error is 

«(*,< + k) - Q h u(x,t) = k[k qi a(x,t) + h q2 b(x,t)} + kO(k qi+1 + h q2+1 ) 

= r -f kO(k qi+1 + h q2+1 ) 
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(5.21) 

(5.22) 



where the leading term is denoted by r. If u is smooth enough, then if one takes two time steps 
with the method Qh, the leading error is 2 r. That is, 

u(se, t + 2k) - Q\u(x, t) = 2r + kO{k^ +l + h^ +1 ) (5.23) 

Let Q 2 h be the same difference scheme as Qh but based on mesh width of 2 h and 2k. Also, assume 
that the order of accuracy in time and space are equal, - q 2 = q. Then 

u(x,t + 2k)-Q 2h u{x,t) = (2k)[(2kya(xJ) + (2hyb{x,t)] + 0(W +2 ) (5.24) 

= 2 9+1 r + 0(h q+2 ) (5.25) 

Since u(x,t + 2k) - Q 2 u(x,i) « 2 r, forming the difference 

Qiu(z,t)-Q 2h u{x.t) = r + 0(/i , +J) ( 5 . 26 ) 

2 ? 2 

gives an estimate of the local truncation error at time t. In other words, the difference between 
the solution on grid (2 h, 2k) and ( h , k) uses to estimate of the local truncation error. 

This procedure has several advantages. First, it is not necessary to know the exact form of 
the truncation error to apply it. The error estimation procedure is independent of the difference 
method. The restriction of this method, that the accuracy in time and space should be the same, 
is not a severe one. Many popular finite difference methods share this property, for example, 
second order methods like Lax Wendroff or MacCormack’s method and Leap Frog, and first order 
method such as upstream differencing. For methods where the accuracy in space and time is not 
the same, a more expensive variant of this procedure is possible. For example, one could estimate 
the spatial and temporal error separately: first keep k constant and take a step based on 2 h 
differences, then keep h constant and take a step with time step 2k. Other variations are possible. 
In the present work the LDS with Richardson extrapolation was employed only for schemes with 
the same spatial and temporal accuracy. The usefulness of this approach applied to discretizations 
without this property should be further investigated. 

The initialization step of the LDS method computes the term Q 2 u — Q 2 hU and uses it as the 
initial value of the correction equation on the next coarser grid. Taking 

(Qlu - Q 2h u) (5.27) 

as the initial value will yield an 0(/z 9+1 ) approximation on the coarse grid. 

The initialization of the extrapolated LDS of degree one is performed directly from the finest 
grid. The extra cost associated with this initialization is compensated by the added accuracy. 

Richardson extrapolation is based on Taylor expansion of the error and is valid only for smooth 
data. For non-smooth solutions, this procedure is incorrect and might lead to an error that grows 
exponentially with the number of visits to the fine grid. Therefore, a great care should be taken 
when considering this method, to ensure that the discretization is dissipative enough to prohibit 
any undesired growth. Nevertheless, since the dissipative schemes employed by the LDS damp the 
oscillatory components (see Sec. 6.2), this technique might give excellent results when properly 
used (see Figure 12). 
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5.3 Work Considerations 


The amount of computational work and the memory requirement in a cycle of an LDS of degree m 
will be evaluated and compared with the corresponding requisites on the finest grid in the cycle. 

The simultaneous time-stepping on the finer grids during initialization constitutes a large 
fraction of the algorithm computational cost and dominates its storage requirements. There is a 
trade off between the storage requirements and the efficiency in the two initialization procedures 
described in Section 5.1. In the present work the computationally efficient simultaneous initial- 
ization was employed since only moderate storage was needed. However, when storage is limited, 
direct initialization might be favored. In this section, only the efficiency of the simultaneous 
scheme is analyzed, but the storage requirements are compared for both methods. 

In order to simplify analysis, it will be assumed that the finest and coarsest grid meshsizes 
satisfy H = 2 ^h\ and that a grid is refined by halving its meshsize. 

The problem is solved in a d-dimensional space, for d = 2,3. Typically, real world problems 
occur in 3-dimensional space. 

5.3.1 Storage Requirements 

The LDS method of degree m employs on the coarsest level m + 1 times as many equations as 
on the finest grid, while on intermediate grids, twice the number of the fine grid equations are 
solved. The number of points in a spatial grid on any level is 2 d larger than in the next coarser 
one. In the simultaneous initialization , see Figure 3, if all grids are allocated simultaneouslv the 
storage requirement is, 

( \ m + 1 L-i 1 \ 

V l + -^r + T,2^iJ S ( 5 . 28 ) 

where S is the storage required for the finest grid. Noting that at all times merely two grids are 

time advanced simultaneously, then a careful management of memory mav reduce this requirement 
to 

f 1 1 ^ + 1\ « 

{ 1 + 2^T + ^dTJ S (5.29) 

These two requisites are equivalent for / < 2 (i.e., for f = 2, or 4) as is the case in the present 
research. Therefore, this possible small reduction of memory usage will not be further elaborated. 

If memory is at premium, storage may be traded for efficiency by using direct initialization , 
see Figure 2. This procedure employs only the finest and coarsest grids with memory requirement 
of 

( 1+ (5.30) 

It should be noted that typically / = 2, d > 2, and due to efficiency considerations the degree 

satisfies m < 3 (see Section 5.3.2); thus, the storage overhead associated with the LDS algorithm 
is fairly small. 

5.3.2 Efficiency 

The computational cost of fine grid time-stepping relative to the cost of obtaining the same 
solution using the LDS algorithm will be evaluated. In this estimate, the cost of the intergrid 
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transfers is neglected, as in many cases it is small relative a fine grid time step. First, the work 
associated with the initialization of an LDS of degree m is computed. Define 

The cost to initialize an LDS of degree to using the simultaneous 
initialization when the spatial dimension is d and H = 2 l h 

An inspection of Figure 3 leads to the following formula for Ii, m (d ), 




= 2m + 


m(m + 1) 

2 dl+1 


/-i j 

+ 3m S o 

k = l z 


(5.32) 


During initialization a time equal to m 2 l fine grid time steps is marched. Denote by N the 
number of coarse grid time steps performed before revisiting the fine grid. The efficiency of an 
LDS cycle is defined as the computational cost of time-stepping on the fine grid relative to the 
work required to obtain the same solution on the coarse grid with an LDS cycle. It is given by 
the formula, 

{N + m )2 33 ', 

h im (d) + N(m+ 1)2-" 

The following table and figure list the LDS efficiency for H = 4h and d — 2, 3, for various values 
of N. 


Revisit 

Efficiency for 2D problems 

Efficiency for 3D problems 


m = 1 

MSEB 

m — 3 


WBM 


5 

6.98 

4.23 

3.24 



4.25 

10 

10.83" 

6.35 

4.67 

16.28 

9.11 

6.63 

15 

13.15 

8.00 




8.83 

20 

15.81 

9.32 

6.75 

27.85 

15.35 

10.86 

25 

17.52 

10.41 






Typically to multilevel methods, this algorithm efficiency increases with the problem dimen- 
sionality. 


6 Fourier Analysis 

Fourier analysis is a major tool for the analysis of numerical approximations to hyperbolic equa- 
tions [17], as well as for analyzing multigrid algorithms [4]. In this section it is employed to obtain 
necessary conditions for algorithm convergence. 

6.1 Properties of the intergrid transfers 

The transfer of the solution between the various grids plays a central role in the LDS algorithm. 
In the initialization stage the solution is first interpolated to the finer grids and after several time 
steps on these grids is restricted back to the coarsest grid. For an LDS of higher degree, this pro- 
cedure is repeated for the correction terms. Inevitably, this process introduces errors. Therefore, 
an appropriate choice of these operators is essential to guarantee the algorithm performances and 
the desired accuracy of the solution. 
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Figure 4: LDS efficiency for 2D and 3D problems when H = Ah. (a) LDS of degree one (b) LDS 
of degree two (c) LDS of degree three. 

The analysis will be first performed for first order operators. This is no limitation since every 
problem may transformed to a first order system by introducing additional variables. Assume 
that A = is fixed on all grids; thus, the temporal error can be expressed in terms of h = 
Aar. For simplicity, the analysis will be performed for one dimensional problems with H = 2h\ 
generalizations to higher dimensions is straightforward. Assume that the spatial discretization is 
of order p. Furthermore, assume that the order of the spatial discretization coincides with the 
order of the full discretization. Thus, for this analysis, the effect of the temporal error on the 
accuracy may be ignored by investigating a semi-discrete system. 

Let the I fa, iff be the interpolation and restriction operators, respectively; and let I ff(0)Jff($) 

be their corresponding symbols. Assume that for smooth components the intergrid operators 
satisfy. 


iff(O) = 1 + c0 Pl + O(0 Pl+1 ) (6.1) 

Iff (6) = 1 + cO Pq + O(0 Pq+1 ) 

where throughout this section c is a generic constant. For the harmonic oscillatory component 
O' = 0 + 7T holds, 


i h H {e') = cd qi + 0(9 qi+1 ) (6.2) 

Iff (O') = cO gq + 0(9 q2+1 ) 

Denote by 1(9) = Iff (6)1 ^(6). Then under the previous assumptions, for smooth components 

1(0) K(l + cO p ')(l + cO p *) = l + c (0 p > +0 p >) + h.o.t (6.3) 
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For the harmonic oscillatory components holds, 



1(9') = c9 c,1+q2 + h.o.t 

(6.4) 

Let L h (6) denote the symbol of L h . For LDS of first degree one wants to compute, 



r(0,A<)= [e Lh W At -e LH W At ]u o (0) 

(6.5) 

Instead one computes, 



[t(0)e Lh{8)At + i(0')e Lh ^ )At - e L ”W At ] u o (0) = t(0, At) + [c (9 Pi + 9 P2 ) e Lh ^ At 

+ c0 <?1+<?2 e £, '‘ (e ' )At ] u 0 (9) + h.o.t 

(6.6) 

For the low frequencies, if L 

is the first derivative, 



8 ( 8 \ 

= * j- + (i) + h °« 

(6.7) 

For the high frequencies, 

= 0 (5) 

(6.8) 

Thus, for this operator, 

t{ 6. At) = c B p+l ~uo(d) + h.o.t 

(6.9) 

In order to obtain the desired accuracy, the following inequalities should hold 



QP1 < QP + 1 ^ 

a 

(6.10) 


gP2 < gp + 1 

h 

(6.11) 


0<7l+<72 < gP+1 ^ 

h 

(6.12) 


In general, for first order hyperbolic equations, ^ = 0(1), thus this term may be neglected and 
one obtains conditions on the order of the intergrid transfers required to guarantee the algorithm 
performances 

P+2 < pi,p 2 (6.13) 

V + 2 < 9i + 92 

It follows from the above bounds that if ^ -*• 0, increasingly higher order interpolations 
will be necessary. Hence, decreasing the time step on a fixed spatial grid will have undesired 
consequences. In practice, one tries to use a time step as large as possible, thus, this observation 
poses no real restrictions. 

If L is an operator of order m, the previous argument implies that, 

f(9, At) = c6 p+m ^u o (0) + h.o.t (6.14) 
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(6.15) 


Therefore, in order to obtain the desired accuracy, the following inequality should hold 

d Pi < 0 p+m — 

h m 

and similar inequalities should hold for P 2 ,(< 7 i + < 72 )- This implies the following conditions, 

p+m+1 < Pl ,p 2 (6.16) 

p + m + 1 < <7i+<72 

provided ££ — 0(1), i.e., it is bounded away from zero. 

For an LDS of degree two, ip is initialized to correct r in the same way that r is initialized to 
correct u LDS . Hence, these conditions are sufficient for LDS of degree two, as well. It follows that 
these conditions are sufficient for an LDS of a general degree m. 

These orders of the intergrid transfer are necessary to maintain the scheme accuracy. In prac- 
tice, one might prefer to employ transfers of order higher than the minimum necessary. Consider, 
for example, a dissipative scheme which strongly damps the oscillatory components. Thus, each 
visit to the fine grid somewhat damps the smoother components of the solution through the inter- 
polation error which transfers some of their energy to the fine grid oscillatory components. High 
order interpolations transfer less energy to these components and therefore might be preferable 
in such circumstances. 

Another undesired property of the intergrid transfer is captured in the following lemma. 
Lemma 2: For an LDS algorithm, if 1(6) + 1(0') = e, with e ^ 1, then initialization procedure 
introduces an 0(1) error in the Fourier component 6. 

Proof : By the assumption, instead of t (0, At) in Eq (6.9) one computes 

+ ( e . /(,)) _ e £»(2„,4<] Ml) = m A() + [ c(jK + 

+ ( € - 1 _ c 0Pi _ c0P 2 )e L^6')At^ Me) + h Q t (g 17) 

The staggered grid transfers, often, satisfy the condition of Lemma 2 for the oscillatory com- 
ponents. However, since the LDS discretizations are, usually, dissipative (see Sec. 6.2), this result 
has little practical consequences. 

6.2 The Symbol of the LDS cycle 

The LDS method employs several operators which should be simultaneously analyzed in order to 
ensure the proper performance of the algorithm. In the sequel, Fourier analysis will be used to 
analyze the cycle of an LDS of degree one for constant coefficient equations in one dimensional 
space with periodic boundary conditions when H — 2h. 

Let Eh, Eh, be the time marching operators on the fine and coarse grids, respectively. Denote 
the intergrid transfers by if? , I fa. The correction term t is initialized by, 

r * = {if! h ~ Eh) u° h (6.18) 

The solution of the LDS of degree one algorithm, for n > 1, is given by, 

U LDS = EjjUfj + (n - 1 )E]f 1 T 1 (6.19) 

= {E H F (n - \)(lf! Ell’fj - Eh)) (6.20) 
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Let E h ,EuJh Jh denote the symbols of the respective operators. Let G(9, u) denote the am- 
plification factor of the Fourier component 9 in a cycle consisting of v consecutive steps on the 
coarse grid before revisiting the fine grid. It is given by 

G{0,v) = II {e h { 20) + * [lF(0)E 2 h (9)ih(6) + i^{9')E 2 h {9')i h H {9') - Eh( 20 )]} ££( 20)11 

The amplification factor of any Fourier component should approximate as closely as possible the 
analytic growth rate of that component. In particular, when the analytic solution does not grow 
in time, one would like to guarantee that no Fourier component is amplified by the LDS cycle, 
i.e., that the amplification factor of the cycle j4(i/) satisfies, 

A(v) = max G(9,v)< 1 (6.21) 

Moreover, it would be desirable if the LDS cycle was time-stable for any underlying discretization 
with this property. ^ 

For scalar equations, it can be easily seen that, if \EhC29 0 )\ = 1, and 

f\29 0 ) = {if?{9o) El(9 0 ) I h H {9 0 ) + 4>o) E 2 h (9' 0 ) I h H {9' 0 ) - E H (29 0 )] u 0 ( 29 0 ) # 0 


then 

lim G(9o,i / ) = oo (6.22) 

z/— * oo 

This is in accordance with the polynomial error bounds derived in Section 3. It should be noted 
that although the bounds were polynomial, whenever condition (6.21) does not hold, the algorithm 
exhibits a growth exponential in the number of visits to the fine grid. Thus, the time-stability 
of the underlying discretization does not necessarily imply time-stability of the LDS scheme. 
Moreover, if this condition is violated for all v < v 0 , (e.g., when the discretization is not dissipative 
enough), then for revisiting index v < i/q the more frequently the fine grid is visited the faster the 
solution blows up. 

Dissipation induces an exponential decay of the solution. Therefore, a dissipative scheme will 
eventually suppress any polynomial growth of the solution, i.e., for such discretizations A(v) < 1 
for large enough v. However, since the LDS algorithm should maintain the fine grid accuracy, 
that grid has to be visited sufficiently often. Therefore, the discretization should have enough 
dissipation to guarantee that condition (6.21) holds for a schedule v prescribed by the accuracy 
requirement. Quite often, artificial dissipation should be added to the coarse grid discretization to 
ensure the cycle is time-stable. This dissipation may be added either during the time marching on 
the coarse grid, or during the initialization stage as well, in which case it should also be introduced 
to the fine grid scheme. 

The additional dissipation typically results in a dissipative cycle. This might limit the method 
applicability for truly long integration times. Hence, one should add just enough dissipation to 
ensure that no Fourier component is amplified by the cycle. This, of course, does not imply that no 
wavenumber may grow during the cycle; to the contrary, imposing such a restrictive requirement 
leads to a dissipative cycle of limited practical interest. 

In general, the length of the cycle v and the amount of artificial dissipation should be deter- 
mined simultaneously to achieve the fine grid accuracy. 
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6.3 Staggered vs Non-staggered grids 

Time dependent problems are often discretized and solved on staggered grids. The dispersion 
error of discretizations on such grids is, typically, significantly smaller than the error in schemes 
of the same order on a non-staggered grid. Hence, one can accurately solve on staggered grids 
for substantially longer integration times. In the sequel the benefits and disadvantages of using 
such grids when employing the LDS method will be briefly discussed. The exposition is rather 
general and is demonstrated through a particular example of the linearized Euler equation, which 
has been investigated in the present research. 

A non-staggered Cartesian grid consists of a set of discrete variables defined at its vertices and 
the system of equations is discretized at these points. In contrast, on a staggered grid the various 
variables are located at different positions in the computational cell and the various equations 
are evaluated at distinct points of the cell. Figure 5 depicts the variables distribution and the 
corresponding computational cells for the linearized Euler equation, 

u t — a • Vu + cp x 

v t = a-Vv + cp y (6.23) 

p t = a ■ Vp + c(u x + v y ) 

where a — (cti,fl 2 ) is ft two dimensional vector. The superiority of staggered grid discretizations 


■ cell for V 


■ cell for U 


□ cell for P 


Figure 5: Staggered grid discretization of the linearized Euler equation 

to those on non-staggered grid stems from the fact that the symbol of a general order mid- 
cell discretization of the first derivative has the form 2 i o i~ ^ — (where a.^ depends 

on the scheme coefficients), while the symbol of a central discretization of this derivative on 
a regular grid is J3/t=i (Pk ftre scheme dependent). The later symbol vanishes at 7 r; 

hence, poorly approximates the continuous symbol in. The phase error of a non-staggered central 
discretization is larger than that of a mid-cell discretization of the same order over a major 
fraction of the spectrum of wavelength representable on the grid. Moreover, in non-dissipative 
full discretizations based on central spatial discretizations the Fourier components in the upper 
half of the spectrum have negative group velocity and move in a direction opposite to the physical 
waves. These components can be spuriously excited by discrete boundary conditions or mesh 
refinements [15, 17]. Therefore, they might deteriorate the computation accuracy even when 
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Figure 6: Refinement of a staggered grid discretization of the linearized Euler equation with zero 
Dirichlet boundary condition in the y direction. 


they do not occur in the physical problem. In contrast, in mid-cell discretizations all Fourier 
components have positive group velocity. A considerable drawback of the staggered grid stems 
from the variables placement at different positions in the computational cell which does not allow 
implementation of characteristic boundary conditions. 

The small phase error of staggered grid discretizations and consequently the significantly less 
artificial dissipation required to suppress the polynomial growth of the LDS error give rise to 
highly efficient algorithms applicable for very long integration time. The efficiency obtained by 
the LDS algorithm for the linear acoustics equation solved on staggered grid was 25 in 2D and 66 
in 3D, (see Figures 15 and 16). 

The staggered grid has two drawbacks associated with the LDS implementation. Both of them 
occur in the model problem investigated here and are associated with the intergrid transfers. The 
first, is that the intergrid transfers satisfy the condition of Lemma 2. Therefore, for some high 
frequencies there is a 0(1) error and enough dissipation should exist to eliminate this error. The 
second problem is that the solution interpolation to the fine grid requires extrapolation in the 
cells nearest to the boundary (see Figure 6). This extrapolation strongly amplifies the oscillatory 
components in the solution and a way had to be found to circumvent this effect for our model 
problem. Despite these deficiencies, for the problems investigated in this work the staggered grid 
supports substantially better results than the regular grid. 


7 Numerical Results 

The LDS method was introduced in the previous sections and its properties were analyzed. The 
numerical examples presented in this section aim at demonstrating the method strength, as well 
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as supporting the claims and analysis performed in previous sections. 

All the examples in this work have spatial dimension two. However, generalization to higher 
dimensional problems is straightforward and was avoided due to the heavy computational cost of 
such simulations. 

In all examples f = 4 with H = ± and h = in Example 6 finer grids were used as well. 

The notation LDS(d = 0,7 = b) used in this section denotes an LDS algorithm of degree a 
which revisits the fine grid after performing b time steps on the coarse grid with the LDS system. 
In all relevant plots, the error is normalized with respect to the solution norm, e.e, the coarse 

u — 11 , where I H U denotes a restriction of the exact solution to 


grid error is normalized by — 
the coarse grid. Similarly, the relative error which measures the error in an approximation to the 
fine grid solution is normalized by the latter solution norm; e.g., the coarse grid relative error is 
normalized by ^ 
grid. 


llJfcjl ^ where u h stands for a restriction of the fine grid solution to the coarse 


7.1 The Advection Equation 

The first set of examples consists of solving the advection equation on the domain [0,1] x [0, 1] 
with periodic boundary conditions. The equation is given by, 

u t - a(x,y)u x - b(x,y)u y = 0 (7.1) 

In case a(x, y),b(x,y) are constant, an explicit solution of this equation with initial data 


u(x,y,0) = u 0 (x,y ) 

(7.2) 

is given by 

u (x, y, t) = U 0 (x + at,y + bt) 

Two instances of this equation are solved, the constant coefficient equation with, 

(7.3) 

a(x,t/) = 1 

b(x,y) = 0.3 

(7.4) 

(7.5) 

and the variable coefficient equation with, 


a(x,y) = 1 + 0.3 sin 2xx 

K x iV) — 0.3 (1 + 0.4 sin 2 ttj/) 

(7.6) 

(7-7) 

In this set of examples, except Example 6, the spatial discretization was second order upwind and 
integration was performed by third order Runge-Kutta. In all these examples £ = K _ 9.3. 

In this set of examples the solution was restricted to the coarse grids by injection, and unless 
otherwise specified, quintic interpolation was employed. 


Example 1. The LDS yields the fine grid solution. Throughout this work, it was claimed 
that the LDS method yields the fine grid solution on the coarse grid. Clearly, these solutions may 
not be identical; rather, one would like to ensure that the error norm in the LDS solution relative 
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to the fine grid solution is similar to the error norm in the fine grid solution relative to the exact 
solution, that is, 

vl ± - _ \\! h u-A\ (78) 

||/f«*|| ~ list'll 

where I h U denotes a restriction of the analytic solution to the fine grid. Figure 7 demonstrates 
that the LDS yields the fine grid solution in the sense defined by (7.8). 

Example 2. The LDS effective integration time. Figure 8 shows two examples of the LDS 
applied to the variable coefficient advection equation, for smooth and oscillatory data. It can be 
readily seen that the effective integration time the discretization can be used (i.e., the integration 
time when the error is small) drastically decreases for oscillatory data. This figure should be born 
in mind as a reference, since in some of the next examples the same equation with these initial 
data are used for integration times significantly longer than the LDS effectiveness time, in order 
to emphasize difference between solutions for various parameters. 

Example 3. Direct and simultaneous initialization are practically indistinguishable. 

In Section 5.1.1 two initialization procedures for the correction terms were presented. These 
methods vary in their computational cost and memory requirements. Figure 9 compares the 
direct and simultaneous initialization procedures. It can be seen that at a time longer than 
the algorithm effectiveness time, the LDS solutions are hardly distinguishable and the observable 
difference between these solutions is significantly smaller than the error in the LDS approximation 
relative to the fine grid solution. Thus, it may be concluded that these procedures have essentially 
the same accuracy. 

Example 4. The necessary and desirable orders of the intergrid transfers. In Sec- 
tion 6.1 the orders of the intergrid transfers required to ensure the algorithm performances were 
analyzed. This analysis implies that for a second order scheme the interpolation should be at 
least cubic for an LDS of general degree d. The injection introduces no error to the high fre- 
quencies, thus is of infinite order. Figure 10 displays an LDS of degree one and two (on the left 
and right, respectively) with linear, cubic and quintic interpolations. It can be seen, for both 
approximations, that linear interpolation results in an approximation worse than the coarse grid 
solution. Cubic interpolation is sufficient to guarantee the LDS performances, however, quintic 
interpolation yields significantly better results. This phenomena might be due to the interpolation 
that transfers some of the smooth components energy to the fine grid high frequencies which are 
strongly damped on the fine grid. This interpolation error is larger for the cubic interpolation. 

Example 5. The efficiency of the LDS of various degrees. An LDS approximation of 
higher degree provides a better approximation to the fine grid solution than a low degree one. 
Hence, it may be used for longer integration time before the fine grid should be revisited. On the 
other hand, for such a scheme the initialization procedure as well as the coarse grid time-stepping 
are more costly. 

Figure 11 addresses the question which approximation is more cost effective, a low or a high de- 
gree one. In this figure LDS algorithms of various degrees with the same efficiencies are compared 
for both smooth and nonsmooth data and for diverse costs. The efficiency of LDS(d = 1,7 = 8), 
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LDS(<f - 2,7 - 20) and LDS(d = 3,7 = 41) is 9.4; while that of LDS(d = 1,7 = 20) and 
LDS(d = 2,7 = 79) is 15.8. In order for an LDS of degree three to achieve efficiency of 15 the 
fine grid should be visited once in 470 coarse grid steps, while an efficiency of 15.8 can not be 
achieved even if the fine grid is visited once in 1000 time steps. This suggests that high efficiency 
can not be achieved with high degree LDS. The efficiency of LDS(d = 3,7 = 350) is only 14.6. 
For short integration times this efficiency can not be achieved with LDS of degree three since not 
enough time steps are marched during the simulation, e.g., see Figure 11 left. 

Figure 11 plots the error in various solutions relative to the fine grid solution. Inspection of 
Figure 11 reveals that for small relative error, e.g. ~ 0.02 (i.e., T < 1 for the smooth data example 
t’h® oscillatory data), then for smooth data the LDS of degree one provides a better 
approximation for both the low and high efficiency scheduling used. For more oscillatory data, 
the degree two LDS yields better results for the same efficiency. This might be explained by the 
fact that LDS of degree one provides a fairly good approximation for the smoother components, 
but not for the more oscillatory ones. The LDS of degree three does not seem an appropriate 
alternative to the lower degree approximations. 

It should be noted that in 3D problems, the performance of LDS of degree two significantly 
improves. Thus, the efficiency of the abovementioned scheduling is : for LDS(d = 1,7 = 8) it 
is 13.63 , for LDS(d = 2,7 = 20) it is 15.34 and for LDS(d = 3,7 = 41) it is 18. Moreover, for 
LDS(d = 1,7 = 20) it is 27.9, and for LDS(d = 2,7 = 79) it is 38.11. Thus, although in 2D 
problems there seems to be little advantage to use second degree LDS rather than first; in 3D, 
this changes drastically. If these schedules yield similar results for a similar equation in 3D, then 
the LDS of degree two is more cost effective for these problems than the first degree algorithm. 

It should be born in mind, that the performances of LDS of various degrees might change de- 
pending on the equations or on the data. Therefore, it is hard to give conclusive recommendations 
which method to prefer; and each problem should be investigated separately. 

Example 6. Richardson Extrapolation. In Section 5.2 the LDS method with Richardson 
extrapolation was introduced. Although this technique is limited in scope and, hence, should 
be used with great care; it might be highly beneficial when it is applicable. Figure 12, provides 
an example when this idea works. It consists of solutions of the constant coefficient advection 
equation discretized with first order upwind forward Euler method, which is first order in time 
and space. The equation is solved on fine grids of 128 and 256 points and corresponding coarse 
grids of 32 and 64 points on which the LDS is solved. In this example the normalized errors are 
computed with respect to the exact solution. For smooth data, mesh refinement of the finest grid 
by a factor of two yields a decrease in error by a factor of 1.875 for the first order scheme versus a 
factor of 5.81 for the LDS with Richardson extrapolation. For nonsmooth data the fine grid error 
is reduced by a factor of 1.71 while for the extrapolated LDS the factor is 4.03. Thus, the error 
reduction for the extrapolated LDS is second order both for smooth and oscillatory data. 

7.2 The Linearized Euler Equation 

The next set of examples involves solution of the linearized Euler equation, given by 

Pt = a(z,y)Px + b(x,y)p y + c(x,y)(u x + v y ) 

u t = a(x,y)u x + b(x,y)Uy + c(x,y)p x 
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(7.9) 



v t = a{x,y)v x + b{x,y)v y + c(x,y)py 

This system was solved for various settings. In all of them interpolation and restriction were of 
sixth order. 

Example 7. Staggered vs. non-staggered grid. Figure 13 demonstrates the superiority of 
the staggered grid discretization over the regular one for the single grid solution and consequently 
for the LDS approximation. It shows solutions of the variable coefficient linearized Euler equation 
on [0, 1] x [0, 1] domain with periodic boundary conditions. The discretization is second order 
upwind for the advection terms and second order central for the terms involving c(z,y); the 
integration was performed with third order Runge-Kutta with ^ = 0.3. In the fine non-staggered 
grid solution one can see an oscillatory component which dominates the solution. This component 
is visible due to a large phase error of the nonsmooth wavelengths on this grid. The error in 
these components is even more visible in the coarse grid and LDS solutions. Clearly, the LDS 
is ineffective for this integration time. However, since the error in the non-staggered fine grid 
solution is already very large, this is not a real drawback. The large dispersive error is also 
the cause for the little resemblance between the staggered and non-staggered solutions for the 
same data and integration time. In contrast, all the staggered grid solutions do not have those 
oscillations, and the LDS provides an excellent approximation to the fine grid solution. 

Example 8. The effect of diminishing CFL. In Section 6.1, it was pointed out that when 
At _ 0, the order of interpolation should increase. Figure 14 demonstrates this phenomenon 
for the constant coefficient linearized Euler equation on [0,1] x [0,1] with periodic boundary 
conditions. Discretization on a staggered grid is second order upwind for the advection term and 
second order central for the terms involving the c factor. Integration was done by low storage 
third order Runge-Kutta. It can be seen that for very short integration time an oscillatory error 
prevails when = 0.01. Increasing by a factor of 10 the CFL as well as the integration time, 
eliminates these oscillations, yielding an excellent approximation to the fine grid solution. 

Example 9. High efficiency LDS on periodic staggered grid. The example in Figure 15, 
is a solution of the constant coefficients acoustics equation, 

Pt ~~ x T Vy 

u t = p x (7-10) 

v t = Py 

on the domain [0, 1] x [0, 1] with periodic boundary conditions discretized on a staggered grid with 
the same discretization as in the previous example, with ^ = 0.3. 

This discretization has only little dissipation through the Runge-Kutta scheme. However, 
since the mid-cell discretization provides an very good approximation to the differential operator, 
the coarse grid operator well approximates the fine grid operator. Hence, only a little artificial 
dissipation should be added. Sixth order artificial dissipation was added to all equations of the 
form c/iA 3 , with e = 0.005. The small dispersive error of the mid-cell discretization enabled both 
very long integration time, as well as very high LDS efficiency of 26 in 2D and 66 in 3D. Note 
that at this stage the LDS error relative to the fine grid is forty times smaller than the relative 
error of the coarse grid. 
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Example 10. High efficiency LDS on non-periodic staggered grid. Figure 16 plots the 
solution of the acoustics equation on [0, 1] x [0, 1] with periodic boundary conditions in the x 
direction for all variables and zero Dirichlet boundary conditions in the y direction for v, 

V(x,0,t) = F(z,l,<) = 0 (7.11) 

except for this difference, all the other parameters are identical to those in the previous example. 

Recall that high order one sided interpolation near the boundary strongly amplifies the oscil- 
latory components. This problem had to be circumvented in this equation for the interpolation in 
the y direction. The observation that led to a resolution of this difficulty is that for this equation 
and these boundary conditions p , u are symmetric in the y direction, and v is asymmetric in this 
direction. These properties were exploited in the design of the intergrid transfers as well as in the 
introduction of the high order artificial dissipation. The p, u variables were svmmetricly extended 
around the boundary in the y direction and the interpolation and dissipation were calculated for 
the extended solution. An assymetric extension for the v variable was similarly defined and used. 

The assymetry v is follows from the dual initial data argument, which asserts that the boundary 
condition (7.11) may be viewed as requiring that a dual solution with the same magnitude but 
opposite sign will constantly hit the boundary from the exterior of the domain (e.g., see [16]). 
The symmetry of p in the y direction follows from 


p y t(x,0,t) = v u (x,0,t) = 0 (7.12) 

Thus, if initially 

Py( x i 0, 0) = 0 (7.13) 

this symmetry is maintained for later time. A similar argument holds for the boundary condition 
at y = 1. For the u variable, 

u yt(x, 0, t) = u t y{x, 0, t) = p xy {x, 0, t ) = p yx (x, 0, t ) = v tx (x, 0, t) = 0 (7.14) 

Again, if the initial solution was symmetric in this direction, symmetry is preserved. These 
assumptions hold for the initial data taken in our examples. 

It can be seen that the efficiency and accuracy of the LDS algorithm were not affected by the 
imposition of non-periodic boundary conditions. 

Example 11. The effect of low order artificial dissipation. The parameters taken in 
the example in Figure 17 are identical to those in Figure 15 except for the use of lower (fourth) 
order of dissipation of the form ehA 2 , with c = 0.02. This choice of c yields the same damping 
of the oscillatory components as the sixth order dissipation used in Example 9. The resulting 
LDS solution does not provide a satisfactory approximation to the fine grid solution even for 
integration time significantly shorter than the one used in Example 9 and more frequent visits 
to the fine grid. It can be concluded that higher order dissipation is indeed essential for the 
algorithm performances, as lower order dissipation damps too strongly the smooth components. 

7.3 The nonlinear Euler equation 

Example 12. The LDS method for the nonlinear Euler equation In Figure 18, the 
nonlinear Euler equation is solved on [0, 1] x [0, 1] domain with periodic boundary conditions. 
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This system is given by, 


Pt 

+ 

p{u x + Uj,) + up x + VPy = 0 



c 2 

Ut 

+ 

(llKj+t) Uy) -\ Px = 0 

Vt 

+ 

c 2 

(u V x + VVy)+ —p y = 0 


(7.15) 


where p = p 1 , c 2 = 7 p' 1 ' -1 and 7 = 1.4. Second order central discretization was used for all terms, 
with third order low storage Runge-Kutta = 0.3). Artificial sixth order dissipation was added 
with e = 0.8. The interpolation and restriction were of sixth order accuracy. 

The LDS efficiency in this example is 17.5 for 2D problems and 32.8 for 3D problems. 


8 Conclusions 

The Large Discretization Step methods for time dependent problems were presented. First, the 
LDS approximation was defined. It consists of a system of lower accuracy discretizations ap- 
proximating a more accurate time dependent discrete operator. Error bounds on this type of 
approximations to linear problems were obtained for the semi-discrete and fully discrete cases. 
These estimates hold for both hyperbolic and parabolic equations. The research reported herein 
aimed at deriving efficient algorithmic implementation of the LDS approximation for hyperbolic 
equations, a type of equations which were not previously amenable to multigrid methods. A 
heuristic argument motivated the design of the LDS algorithm for nonlinear problems, as well. 

The LDS methods enables to obtain the fine grid accuracy on a coarse grid by adding correction 
terms to the coarse grid equations, initializing them using the fine grid and solving a system of 
equations for these terms. The accuracy of the correction terms deteriorates at a rate determined 
by the coarse grid discretization. However, since their norm is significantly smaller than the 
solution norm, they may be effectively used for many coarse grid time steps. Thereafter, the fine 
grid should be revisited to compute new initial data for them. Fourier analysis was employed to 
analyzed different aspects of the algorithm; in particular, to obtain conditions on the necessary 
orders of the intergrid transfers. 

The resulting algorithm has a typical efficiency of 16 for 2D problems and 28 for 3D equations. 
This efficiency was achieved for linear problems with periodic and Dirichlet boundary conditions 
and the for the nonlinear Euler equation with periodic boundary conditions. A particularly good 
discretization of a linear equation yielded efficiency of 25 in 2D and 66 in 3D problem. 

The results presented in this work are very promising as for the potential of the proposed 
approach to tackle even more complex problems. Still, a lot of research should be carried out 
to better understand the method abilities and limitations. Several aspects of the LDS algorithm 
should be further investigated. The frequency the fine grid is visited is important both to the 
algorithm efficiency and to ensure the accuracy of the resulting solution. It would be highly ben- 
eficial if there was a systematic, preferably adaptive, way to determine when the fine grid should 
be visited. The dissipativity of the LDS cycle is essential to damp the polynomially growing error. 
For dissipative discretizations, this requirement typically does not pose any problems. Quite of- 
ten, though, some artificial dissipation should be added to guarantee the algorithm performances. 
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A general way is required to determine the right amount of artificial dissipation which will sup- 
press the error growth without affecting the algorithm accuracy. Boundary conditions treatment 
within the LDS does not seem to pose particular problems. Nevertheless, -a large variety of such 
conditions should be implemented and tested. In particular, a general approach should be found 
to reduce the high interpolation error near the boundaries. 

The novel character of these methods opens many research directions. The application of the 
LDS to systems of conservations laws and shocks calculations should be investigated. Another 
interesting direction is to investigate possible generalizations of this method. In the present 
research, one solves for the correction terms of a linear problem the same equation as the solution. 
This approximation does not always yield the desired results, e.g., consider the discretization of 
the wave equation discussed in Section 3. In a more general setting, there might be several 
correction terms each satisfying a different equation. Such a generalization might significantly 
reduce the simplicity of the present approach; however, it could be applicable to a broader class 
of equations. 

It is expected that the incorporation of the LDS ideas into parabolic solvers would significantly 
improve their performances, as well. This is another promising research direction. 
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Figure 7: The advection equation, u t + u x + 0.3 u y = 0, with initial data u 0 = e -20(* 2 +s , 2 ) 
Left : (a) fine grid solution (b) the LDS(d = 1, 7 = 20) solution (c) coarse grid solution 
Right : Relative errors, (a) \\I h u - t» A ||/||J*t t ||, (b) \\u h - u» DS ||/||t* A ||, (c) ||« A - t t »||/|| tt *||. 



Figure 8: Solutions of u t + (1 + 0.3 sin(2jr*))u r + 0.3 (1 + 0.4 cos(2iry)) Uj/ = 0. Left : Initial data 
is u 0 = e <* +* ). Right : Initial data is u 0 = e - 50 (^+v 2 ). I n both figures : (a) the fine grid 
solution (b) the LDS(d = 1,7 = 20) solution (c) the coarse grid solution. 







T = 3.2 



T = 3.2 



Figure 9: The advection equation, u t + (1 + 0.3 sin(27rx))u x + 0.3(1 + 0.4 cos(2ir y))u y = 0, with 
Uo = e -20(x 2 +y 2 ) L e f t : ( a ) the fine grid solution (b) the LDS(d = 2,7 = 20) solution when 
employing direct initialization, (c) the coarse grid solution (d) the LDS(d = 2,7 = 20) solution 
when employing simultaneous initialization. Right : Relative errors, (b) relative error of the 
LDS(d = 2,7 = 20) solution when employing direct initialization, (c) relative error of the coarse 
grid solution (d) relative error of the LDS(d = 2,7 = 20) solution when employing simultaneous 
initialization. 
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Figure 10: Solution of u<+(l+0.3 sin(27ra:))u r +0.3 (1+0.4 cos(27ry))u y = 0,with u 0 = e~ 20 i^+v 2 ). 
Left : (a) the fine grid solution (b) the LDS(d =1,7 = 20) solution when using quintic interpola- 
tion (c) coarse grid solution (d) the LDS(d = 1,7 = 20) solution when using cubic interpolation 
(e) the LDS(d = 1,7 = 20) solution when using linear interpolation. Right : (a) the fine grid solu- 
tion (b) the LDS(d = 2,7 = 20) solution when using quintic interpolation (c) coarse grid solution 
(d) the LDS(d = 2,7 = 20) solution when using cubic interpolation (e) the LDS(d = 2,7 = 20) 
solution when using linear interpolation. 
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Figure 11: Relative error in the solution of u t + (l + 0.3 sin(27rar))tt I + 0.3 (1 + 0.4 cos(2^y))u y — 0, 
with Left: Initial data is u 0 = € -20(* 2 +y s ) Right : Initial data is u 0 = e ~ 50(x2+y2) . In both 
figures : (a) the LDS(d = 1,7 = 8) solution (b) the LDS(d = 2,7 = 20) solution (c) the 
LDS(d = 3,7 = 41) solution (d) the LDS(d = 1,7 = 20) solution (e) the LDS(d = 2,7 = 79) 
solution (f) the LDS(d = 3,7 = 350) solution. 
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Figure 12: Errors in solutions of u t + u x + 0.3 u y = 0. Left : Initial data is u 0 = e -25(r= + y 2 ) 
Right : Initial data is u 0 = e -60(x 2 + y 2 ) In both figures . ^ error in so ] ution on a 12g x 12 g 

points grid, (b) error in LDS(d = 1,7 = 20) solution on a 32 x 32 points grid with Richardson 
extrapolation, (c) error in solution on a 256 x 256 points grid, (d) error in LDSfd = 1,7 = 20) 
solution on a 64 x 64 points grid with Richardson extrapolation. 


47 





Figure 13: Solution of the variable coefficient linearized Euler equation, with coefficients : 
a(x v) = 0.6(1 + 0.2 cos27r:c), b(x,y) = 0.4(1 + 0 . 4 cos 2 ttj /), c(x,y) = 1 + 0.3sm27T.r. Ini- 
tialdata is On = e - 30 ^ 2+y2) , tto = 0, v 0 = 0. Left : non-staggered grid. Right : staggered grid. 
In both figures : (a) the fine grid solution (b) the LDS(d = 1.7 = 20) solution (c) the coarse grid 

solution. 



Figure 14: Solution of the constant coefficients linearized Euler equation, with coefficients : 
a = 0.4, b = 0.5, c = 1. Initial data is p 0 = e " 45(z2+y2) , u 0 = 0, v 0 = 0. Left : ^ = 0.01. 
Right - A = 0.1. In both figures : (a) the fine grid solution (b) the LDS(d = 1,7 = 2 0 ) solution 

(c) the coarse grid solution. 







Figure 15: The acoustics equation on periodic domain, with initial data p 0 = e~ 35 ^ x2+y2 \ uq = 0 
v 0 = 0. In the plots of the solutions : (a) the fine grid solution (b) the LDS(d = 1, 7 = 80) solution 
when adding sixth order dissipation (c) the coarse grid solution. In the plot of the relative error • 

(b) the relative error in the LDS(d = 1, 7 = 80) solution (c) the relative error in the coarse grid 
solution. ° 







Figure 16 - Solution of the acoustics equation with po - e 35 ^ +y * , u 0 — 0, u 0 - 0 on non- 
periodic domain, (a) the fine grid solution (b) the LDS(d = 1,7 = 80) solution when adding sixth 
order dissipation (c) the coarse grid solution. 
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